#############################
create_list_expected_profits_2018<- function(TM2,df){
  #df are the product chosen (i.e borrower info), TM is the data in the right format
  #predict the product share, amount borrowed and default
  
  ###############################
  
  predict_logit<-predict(model2,newdata=TM2) #matrix with borrower id as rows and predicted MS as colums (IDTBP)
  predict_logit<-as.data.frame(predict_logit)
  coef_logit<-as.data.frame(t(coef(model2)))
  rate_coef<-coef_logit$initial_gross_interest_rate + df$gross_income*coef_logit$`initial_gross_interest_rate:gross_income`
  
  coef_default<-as.data.frame(t(coef(modelDefault)))
  rate_default_coef<-coef_default$initial_gross_interest_rate
  
  # Do those predictions for each product-bank
  product_data<-unique(df[,c("IDTBP","Fee","Max_LTV","initial_gross_interest_rate","Teaser_rate_years","BANKID")])
  
  #initialize my.list.expected.profit 
  ######
  my.list.expected.profit<-list()
  borrowerID<-0
  BANKID<-0
  IDTBP<-0
  predict_Loan<-0
  predict_Default<-0
  share<-0
  initial_gross_interest_rate<-0
  Max_LTV<-0
  Fee<-0
  Teaser_rate_years<-0
  for(b in sort(unique(product_data$BANKID)) ){
    my.list.expected.profit[[b]]<-data.frame(borrowerID,BANKID,IDTBP,predict_Loan,predict_Default,share,initial_gross_interest_rate,Max_LTV,Fee,Teaser_rate_years)
  }
  
  
  #######
  
  
  i<-0 #used to get the right IDTBP from the predict_logit
  for(idtbp in sort(unique(product_data$IDTBP))){
    i<-i+1
    df_bis<-df
    df_bis$Fee<-product_data$Fee[product_data$IDTBP==idtbp]
    df_bis$Max_LTV<-product_data$Max_LTV[product_data$IDTBP==idtbp]
    df_bis$initial_gross_interest_rate<-product_data$initial_gross_interest_rate[product_data$IDTBP==idtbp]
    df_bis$Teaser_rate_years<-product_data$Teaser_rate_years[product_data$IDTBP==idtbp]
    df_bis$BANKID<-product_data$BANKID[product_data$IDTBP==idtbp]
    df_bis$IDTBP<-idtbp
    
    predict_Loan<-exp(predict(modelLoan,newdata=df_bis))
    predict_Default<-predict(modelDefault,newdata=df_bis)
    
    initial_gross_interest_rate<-df_bis$initial_gross_interest_rate
    Max_LTV<-df_bis$Max_LTV
    Fee<-df_bis$Fee
    Teaser_rate_years<-df_bis$Teaser_rate_years
    
    borrowerID<-df_bis$borrowerID
    BANKID<-df_bis$BANKID
    IDTBP<-df_bis$IDTBP
    share<-predict_logit[,i]
    
    tmp<-as.data.frame(rbind(my.list.expected.profit[[BANKID[1]]],data.frame(borrowerID,BANKID,IDTBP,predict_Loan,predict_Default,share,initial_gross_interest_rate,Max_LTV,Fee,Teaser_rate_years)))
    my.list.expected.profit[[BANKID[1]]]<-tmp
  }
  
  
  my.list.expected.profits<-list()
  for (b in 1:max(df$BANKID)) {
    my.list.expected.profits[[b]]<-my.list.expected.profit[[b]][2:(nrow(my.list.expected.profit[[b]])-1),]
    # my.list.expected.profits[[b]]<-my.list.expected.profits[[b]][!is.na(my.list.expected.profits[[b]]$IDTBP),]
  }
  
  return(list(my.list.expected.profits,rate_coef,rate_default_coef))
}
#############################

#############################
create_Matrices_2018 <- function(j,my.list.expected.profit.alpha,alpha_d,mean_alpha){ # Work, to add, the maturity of the period (or the fix rate duration)
  #my.list.expected.profit.alpha<-list_a 
  # alpha_d: coefficient in front of rate in the default
  # mean_alpha<-rate_coef: coefficient in front of rate in the logit (different for each borrower)
  sigma<-1 # scaling parameter (not needed here)
  A_C_jt<-list()
  A_C_jt2<-list()
  
  B_C_jt<-list()
  B_C_jt2<-list()
  B_C_jt3<-list()
  
  R_list<-list() # drop the  NA
  R_couner_list<-list() # contains NA
  Ones_list<-list() # drop the  NA
  
  Max_LTV_list<-list()
  Fee_list<-list()
  TEASER_list<-list()
  
  #####################################################################
  # for (counter in 1:6) {
  counter<-4
  # start.time<- Sys.time()
  # print(start.time)
  
  # browser()
  
  nb_prod<-length(unique(my.list.expected.profit.alpha[[1]]$IDTBP)) 
  # print(counter)
  print(nb_prod)
  
  if (is.na(nb_prod) ){ # if the list is empty (i.e. that bank does not offer product of that type counter type)
    # print(counter)
    A_C_jt[[counter]]<-NULL
    A_C_jt2[[counter]]<-NULL
    
    B_C_jt[[counter]]<-NULL
    B_C_jt2[[counter]]<-NULL
    B_C_jt3[[counter]]<-NULL
    R_list[[counter]]<-NULL
    Ones_list[[counter]]<-NULL
    R_couner_list[[counter]]<-NA
    Max_LTV_list[[counter]]<-NULL
    Fee_list[[counter]]<-NULL
    TEASER_list[[counter]]<-NULL
    
  }else{
    #loop for each alpha 
    
    # Create a matrix A (alpha) [contract number j][ derivative with respcted to r_j]
    # A_xy_alpha<-list(list())
    A_xy_alpha<-matrix(NA,nb_prod,nb_prod)
    A_xy_alpha2<-matrix(NA,nb_prod,nb_prod)
    B_xy_alpha<-matrix(NA,nb_prod,nb_prod)
    B_xy_alpha2<-matrix(NA,nb_prod,nb_prod)
    B_xy_alpha3<-matrix(NA,nb_prod,nb_prod)
    R<-NULL
    Max_LTV<-NULL
    Fee<-NULL
    TEASER<-NULL
    
    nb_alpha<-1 ####### assume just 3 draws for simplicity
    for (alpha in 1:nb_alpha) {
      # print(alpha)
      dff<-my.list.expected.profit.alpha[[alpha]] 
      
      # loop for each product
      idmatprod1<-0
      for (idxprod in sort(unique(my.list.expected.profit.alpha[[1]]$IDTBP)) ) {
        idmatprod1<-idmatprod1+1
        # print(idxprod)
        subdf<-dff[dff$IDTBP==idxprod,] #should contain only 3640 observations
        # ixd_S<- which(names(dff)%in%c(paste0("s_alpha_",idxprod)))
        # ixd_D<- which(names(dff)%in%c(paste0("D_alpha_",idxprod)))
        # ixd_L<- which(names(dff)%in%c(paste0("L_alpha_",idxprod)))
        # ixd_r<- which(names(dff)%in%c(paste0("initial_gross_interest_rate_",idxprod)))
        # idx_Max_LTV<-which(names(dff)%in%c(paste0("Max_LTV_",idxprod)))
        
        # Mat<-25 ######################################################to add, the maturity of the period (or the fix rate duration)
        
        Mat<-25 ######################################################to add, the maturity of the period (or the fix rate duration)
        
        
        S_alpha<-subdf$share
        D_alpha<-subdf$predict_Default
        L_alpha<-subdf$predict_Loan
        r<-subdf$initial_gross_interest_rate
        loantoValue<-subdf$Max_LTV
        Fee<-subdf$Fee
        teaser<-subdf$Teaser_rate_years
        
        
        rate<-r
        
        # colnames(rate)[1] <- "rate"
        # colnames(loantoValue)[1]<- "loantoValue"
        # browser()
        
        
        
        
        if (alpha==1) {
          # browser()
          R<-c(R,rate[1]) #a<-rstore the rates
          Max_LTV<-c(Max_LTV,loantoValue[1]) #a<-rstore the Max_LTV
          Fee<-c(Fee,Fee[1])
          TEASER<-c(TEASER,teaser[1])
          
          idmatprod2<-0
          for (idxprod2 in sort(unique(my.list.expected.profit.alpha[[1]]$IDTBP))) {
            idmatprod2<-idmatprod2+1
            M_S_alpha<-rep(0, nb_prod) 
            
            if (idmatprod1==idmatprod2) {
              # browser()
              # sum( (mean_alpha+ B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha*(1+sigma)*Mat) # sum for all borrowers 
              
              #the diagonal elements
              A_xy_alpha[idmatprod1,idmatprod2]<- sum( (mean_alpha+ B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha*(1+sigma)*Mat) # sum for all borrowers 
              A_xy_alpha2[idmatprod1,idmatprod2]<- sum( sigma*(mean_alpha+ B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha*S_alpha*Mat) # sum for all borrowers 
              
              B_xy_alpha[idmatprod1,idmatprod2]<- sum( B_i_r_weight[alpha]*S_alpha*L_alpha *Mat*(  (mean_alpha+ B_i_r[alpha])* (  (1-D_alpha)* (1+sigma) )+ (   - alpha_d)  )  )   # sum for all borrowers 
              B_xy_alpha2[idmatprod1,idmatprod2]<- sum( B_i_r_weight[alpha]*S_alpha*L_alpha *Mat*(  sigma*(mean_alpha+ B_i_r[alpha])* (  (1-D_alpha)* ((S_alpha)) ) )  )   # sum for all borrowers 
              B_xy_alpha3[idmatprod1,idmatprod2]<- sum( B_i_r_weight[alpha]*S_alpha*L_alpha *Mat*(   ((1-D_alpha))   )  )   # Q 
              
              
              # B_i_off[[idxprod2]]<-0
            }else{
              subdf2<-dff[dff$IDTBP==idxprod2,]
              S_alpha2<-subdf2$share
              #the off diagonal elements (derivative with respect to r_j jdifferent than idx prod)
              A_xy_alpha[idmatprod1,idmatprod2]<- 0
              A_xy_alpha2[idmatprod1,idmatprod2]<- sum( sigma*(mean_alpha+ B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha* S_alpha2*Mat) #get market share
              
              B_xy_alpha[idmatprod1,idmatprod2]<- 0
              B_xy_alpha2[idmatprod1,idmatprod2]<-sum(sigma*(mean_alpha+ B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha*Mat*(1-D_alpha)* S_alpha2)
              B_xy_alpha3[idmatprod1,idmatprod2]<-0
            }
            
          }
          
          
        }else{
          
          idm2<-0
          for (idxprod2 in sort(unique(my.list.expected.profit.alpha[[1]]$IDTBP))) {
            idmatprod2<-idm2+1
            M_S_alpha<-rep(0, nb_prod) 
            
            if (idmatprod2==idmatprod1) {
              # browser()
              A_xy_alpha[idmatprod1,idmatprod2]<-A_xy_alpha[idmatprod1,idmatprod2] +  sum( (mean_alpha+ B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha*(1+sigma)*Mat) # sum for all borrowers 
              A_xy_alpha2[idmatprod1,idmatprod2]<-A_xy_alpha2[idmatprod1,idmatprod2]+ sum( sigma*(mean_alpha+ B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha*S_alpha*Mat) # sum for all borrowers 
              
              B_xy_alpha[idmatprod1,idmatprod2]<-B_xy_alpha[idmatprod1,idmatprod2]+ sum( B_i_r_weight[alpha]*S_alpha*L_alpha *Mat*(  (mean_alpha+ B_i_r[alpha])* (  (1-D_alpha)* (1+sigma) )+ (   - alpha_d)  )  )   # sum for all borrowers 
              B_xy_alpha2[idmatprod1,idmatprod2]<-B_xy_alpha2[idmatprod1,idmatprod2]+  sum( B_i_r_weight[alpha]*S_alpha*L_alpha *Mat*(  sigma*(mean_alpha+ B_i_r[alpha])* (  (1-D_alpha)* ((S_alpha)) ) )  )   # sum for all borrowers 
              B_xy_alpha3[idmatprod1,idmatprod2]<- B_xy_alpha3[idmatprod1,idmatprod2]+  sum( B_i_r_weight[alpha]*S_alpha*L_alpha *Mat*(   ((1-D_alpha))   )  )   # sum for all borrowers 
              
              
            }else{
              #the off diagonal elements (derivative with respect to r_j jdifferent than idx prod)
              subdf2<-dff[dff$IDTBP==idxprod2,]
              S_alpha2<-subdf2$share
              
              A_xy_alpha2[idmatprod1,idmatprod2]<-  A_xy_alpha2[idmatprod1,idmatprod2] + sum(sigma*(mean_alpha+B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha* S_alpha2*Mat) #get market share
              
              B_xy_alpha2[idmatprod1,idmatprod2]<-B_xy_alpha2[idmatprod1,idmatprod2] + sum(sigma*(mean_alpha+ B_i_r[alpha])*B_i_r_weight[alpha]*S_alpha*L_alpha*Mat*(1-D_alpha)* S_alpha2)
              
            }
          }
        }
      }
    }
    
    # print(counter)
    A_C_jt[[counter]]<-A_xy_alpha
    A_C_jt2[[counter]]<-A_xy_alpha2
    
    
    B_C_jt[[counter]]<-B_xy_alpha
    B_C_jt2[[counter]]<-B_xy_alpha2
    B_C_jt3[[counter]]<-B_xy_alpha3
    
    R_list[[counter]]<-R
    R_couner_list[[counter]]<-R
    
    Ones_list[[counter]]<-rep(1,length(R))
    
    Max_LTV_list[[counter]]<-Max_LTV
    Fee_list[[counter]]<-Fee
    TEASER_list[[counter]]<-TEASER
  }
  # end.time <- Sys.time()
  # print(end.time)
  # time.taken <- end.time - start.time
  # print(time.taken)
  # }
  
  # DO the matrix for the B part as well
  
  # drop A_C_jt when the value is NULL
  
  if (length(A_C_jt)>=1) {
    # browser()
    A_C_jt<-A_C_jt[!sapply(A_C_jt, is.null)]
    A_C_jt2<-A_C_jt2[!sapply(A_C_jt2, is.null)]
    
    B_C_jt<-B_C_jt[!sapply(B_C_jt, is.null)]
    B_C_jt2<-B_C_jt2[!sapply(B_C_jt2, is.null)]
    B_C_jt3<-B_C_jt3[!sapply(B_C_jt3, is.null)]
    
    R_list<-R_list[!sapply(R_list, is.null)]
    Max_LTV_list<-Max_LTV_list[!sapply(Max_LTV_list, is.null)]
    Fee_list<-Fee_list[!sapply(Fee_list, is.null)]
    TEASER_list<-TEASER_list[!sapply(TEASER_list, is.null)]
    Ones_list<- Ones_list[!sapply(Ones_list, is.null)]
    
    ##
    A_jt<-.bdiag(A_C_jt) # tilde Lambda
    A_jt2<-.bdiag(A_C_jt2)# tilde Gamma
    
    B_jt<-.bdiag(B_C_jt) #  Lambda
    B_jt2<-.bdiag(B_C_jt2)#  Gamma
    B_jt3<-.bdiag(B_C_jt3) # Q
    
    RR<-.bdiag(R_list) # R
    Max_LTVMax_LTV<- Max_LTV_list #.bdiag() # Max_LTV
    FeeFee<- Fee_list #.bdiag() # Max_LTV
    TEASERTEASER<- TEASER_list #.bdiag() # Max_LTV
    ONES<- .bdiag(Ones_list) # 1
    
    # print(R_couner_list)
    
    # return(solve(A_jt, B_jt))
    return(list(A_jt,A_jt2, B_jt,B_jt2,B_jt3,RR,R_couner_list,Max_LTVMax_LTV,ONES,FeeFee,TEASERTEASER))
    
  }else {
    
    return(NA)
    
  }
  ### problem if singular matrix, #may need to return the characteristics of the product
  
} #############################

#############################
create_list_TM_df_new_prod <- function(TM2,merged_menu,subdt){
  #Add the new product in the borrower choice set
  
  m<-merged_menu[merged_menu$offered_products==0,]
  
  TM2NEW.l<-list()
  dfNEW.l<-list()
  # nrow(m)
  
  # browser()
  for (i in 1:80) { # normally put 1:nrow(m), put too big atm
    subdtNew<-subdt[subdt$choice==1,] # dataframe to replace
    subdtNew$IDTBP<-m$IDTBP[i]
    subdtNew$Teaser_rate_years <-m$Teaser_rate_years[i]
    subdtNew$Fee<-m$round_fees[i]
    subdtNew$Max_LTV<-m$round_LTV[i]
    subdtNew$initial_gross_interest_rate <-m$initial_gross_interest_rate[i]
    subdtNew$BANKID  <-m$BANKID[i]
    subdtNew$choicevar<-m$IDTBP[i]
    subdtNew$choice<-0
    subdtNew<-rbind(subdt,subdtNew)
    TM2NEW.l[[i]]<-mlogit.data(subdtNew, choice = "choice", shape = "long", chid.var = "borrowerID", alt.var = "IDTBP")
    
    #chose one borrower and swap his choice to the new product
    dd<-as.data.frame(TM2NEW.l[[i]][TM2NEW.l[[i]]$choice==1,])
    dd$Teaser_rate_years[1] <-m$Teaser_rate_years[i] #replace the first  row by a new product
    dd$Fee[1]<-m$round_fees[i]
    dd$Max_LTV[1]<-m$round_LTV[i]
    dd$initial_gross_interest_rate[1] <-m$initial_gross_interest_rate[i]
    dd$BANKID[1]  <-m$BANKID[i]
    dd$IDTBP[1]  <-m$IDTBP[i]
    dfNEW.l[[i]]<-dd
  }
  return(list(TM2NEW.l,dfNEW.l))
}

#############################

#############################
create_list_TM_df_no_high_Max_LTV <- function(TM2,merged_menu,subdt){
  #create the IDTBP of the new product
  # IDTBPNEW<-max(as.integer(as.vector(unique(TM2$IDTBP))))+1
  
  # m<-merged_menu[merged_menu$offered_products==0,]
  # 
  TM2NEW.l<-NULL
  dfNEW.l<-NULL
  # nrow(m)
  
  # browser()
  i<-1
  subdtNew<-subdt # dataframe to replace
  # dataset of low Max_LTV products
  subdtlow_Max_LTV<-subdtNew[subdtNew$Max_LTV<92,]
  
  #the next line will create duplicate products
  subdtNew$duplicate<-0
  subdtNew[subdtNew$Max_LTV>=92&subdtNew$choice=="TRUE",]$duplicate<-1
  #replace the choice of all borrowers choosing high Max_LTV to the first contract of subdtlow_Max_LTV
  subdtNew[subdtNew$Max_LTV>=92&subdtNew$choice=="TRUE",c("IDTBP","Teaser_rate_years","Fee","Max_LTV","initial_gross_interest_rate","BANKID")]<-subdtlow_Max_LTV[1,c("IDTBP","Teaser_rate_years","Fee","Max_LTV","initial_gross_interest_rate","BANKID")]
  #get rid of the duplicate
  subdtNew<-subdtNew[!(subdtNew$borrowerID %in% subdtNew$borrowerID[subdtNew$duplicate==1] & subdtNew$IDTBP==subdtlow_Max_LTV[1,c("IDTBP")] & subdtNew$choice=="FALSE"),]
  # keep only the menus with low Max_LTV
  subdtNew<-subdtNew[subdtNew$Max_LTV<92,]
  
  # browser()
  TM2NEW.l<-mlogit.data(subdtNew, choice = "choice", shape = "long", chid.var = "borrowerID", alt.var = "IDTBP")
  
  #chose one borrower and swap his choice to the new product
  dd<-as.data.frame(TM2NEW.l[TM2NEW.l$choice=="TRUE",])
  dfNEW.l<-dd
  # dd$Teaser_rate_years[1] <-m$Teaser_rate_years[i] #replace the first  row by a new product
  # dd$Fee[1]<-m$round_Fee[i]
  # dd$Max_LTV[1]<-m$round_Max_LTV[i]
  # dd$initial_gross_interest_rate[1] <-m$initial_gross_interest_rate[i]
  # dd$BANKID[1]  <-m$BANKID[i]
  # dd$IDTBP[1]  <-IDTBPNEW
  # dfNEW.l[[i]]<-dd
  
  return(list(TM2NEW.l,dfNEW.l))
}
#############################

#############################
create_list_TM_df_no_high_Max_LTV_add_product <- function(TM2,prod_add,subdt,nb_prod){
  #create the IDTBP of the new product
  IDTBPNEW<-max(as.integer(as.vector(unique(TM2$IDTBP))))
  
  # subdtb
  
  #product to be added 
  prod_add$IDTBP<-IDTBPNEW+1
  
  
  TM2NEW.l<-list()
  dfNEW.l<-list()
  # nrow(m)
  
  subdtNew<-subdt # dataframe to replace
  # dataset of low Max_LTV products
  
  
  
  #add the new product 
  # browser()
  subdtN<-NULL
  # normally put 1:nrow(m), put too big atm
  # browser()
  IDTBPNEW<-IDTBPNEW+1
  subdtNew<-subdt[subdt$choice=="TRUE",] # dataframe to replace
  subdtNew$IDTBP<-rep(as.numeric(IDTBPNEW),nrow(subdtNew))
  subdtNew$Teaser_rate_years <-rep(prod_add$Teaser_rate_years,nrow(subdtNew))
  subdtNew$Fee<-rep(prod_add$Fee,nrow(subdtNew))
  subdtNew$Max_LTV<-rep(prod_add$Max_LTV,nrow(subdtNew))
  subdtNew$initial_gross_interest_rate <-rep(prod_add$initial_gross_interest_rate,nrow(subdtNew))
  subdtNew$BANKID  <-rep(prod_add$BANKID,nrow(subdtNew))
  subdtNew$choice<-rep("FALSE",nrow(subdtNew))
  # browser()
  #create the dataset with the new product
  subdtNew<-rbind(subdt,subdtNew)
  
  # browser()
  #chose  borrower nb_prod and swap his choice to the new products (otherwise logit estimation bug)
  a<-subdtNew[subdtNew$choice=="TRUE",]
  borrowerID<-a[nb_prod,]$borrowerID
  subdtNew[subdtNew$choice=="TRUE"&subdtNew$borrowerID==borrowerID,]$choice<-"FALSE"
  subdtNew[subdtNew$IDTBP==as.numeric(IDTBPNEW)&subdtNew$borrowerID==borrowerID,]$choice<-"TRUE"
  
  
  TM2NEW.l<-mlogit.data(subdtNew, choice = "choice", shape = "long", chid.var = "borrowerID", alt.var = "IDTBP")
  
  dfNEW.l<-as.data.frame(TM2NEW.l[TM2NEW.l$choice=="TRUE",])
  
  
  # browser()
  #chose borrowers that chose the same contracts and swap his choice to the new products
  # dd<-as.data.frame(TM2NEW.l[TM2NEW.l$choice=="TRUE",])
  # a<-summary(as.factor(dd$IDTBP))
  # ID_ML<-as.numeric(row.names(as.matrix(a[1]))) # get the numeric value of the most likely occurrence f IDTBP
  # 
  # IDTBPNEW<-max(as.integer(as.vector(unique(TM2$IDTBP))))
  # for (i in 1:nrow(m)) {
  #   IDTBPNEW<-IDTBPNEW+1
  #   dd[dd$IDTBP==ID_ML,]$Teaser_rate_years[i] <-m$Teaser_rate_years[i]
  #   dd[dd$IDTBP==ID_ML,]$Fee[i]<-m$round_Fee[i]
  #   dd[dd$IDTBP==ID_ML,]$Max_LTV[i]<-m$round_Max_LTV[i]
  #   dd[dd$IDTBP==ID_ML,]$initial_gross_interest_rate[i]<-m$initial_gross_interest_rate[i]
  #   dd[dd$IDTBP==ID_ML,]$BANKID[i] <-m$BANKID[i]
  #   dd[dd$IDTBP==ID_ML,]$IDTBP[i] <-IDTBPNEW
  # }
  # 
  # dfNEW.l<-dd
  
  return(list(TM2NEW.l,dfNEW.l))
}
#############################

#############################
create_list_TM_df_no_high_Max_LTV_add_teaser <- function(TM2,merged_menu,subdt){
  #create the IDTBP of the new product
  IDTBPNEW<-max(as.integer(as.vector(unique(TM2$IDTBP))))
  
  #may need to change this
  m<-merged_menu[merged_menu$offered_products==0 & merged_menu$round_Max_LTV<90 &merged_menu$Teaser_rate_years==7 ,]
  
  TM2NEW.l<-list()
  dfNEW.l<-list()
  # nrow(m)
  
  subdtNew<-subdt # dataframe to replace
  # dataset of low Max_LTV products
  subdtlow_Max_LTV<-subdtNew[subdtNew$Max_LTV<90,]
  
  #the next line will create duplicate products
  subdtNew$duplicate<-0
  subdtNew[subdtNew$Max_LTV>=90&subdtNew$choice=="TRUE",]$duplicate<-1
  #replace the choice of all borrowers choosing high Max_LTV to the first contract of subdtlow_Max_LTV
  subdtNew[subdtNew$Max_LTV>=90&subdtNew$choice=="TRUE",c("IDTBP","Teaser_rate_years","Fee","Max_LTV","initial_gross_interest_rate","BANKID")]<-subdtlow_Max_LTV[1,c("IDTBP","Teaser_rate_years","Fee","Max_LTV","initial_gross_interest_rate","BANKID")]
  #get rid of the duplicate
  subdtNew<-subdtNew[!(subdtNew$borrowerID %in% subdtNew$borrowerID[subdtNew$duplicate==1] & subdtNew$IDTBP==subdtlow_Max_LTV[1,c("IDTBP")] & subdtNew$choice=="FALSE"),]
  # keep only the menus with low Max_LTV
  subdtb<-subdtNew[subdtNew$Max_LTV<90,]
  
  #add the new products 
  # browser()
  subdtN<-NULL
  for (i in 1:nrow(m)) { # normally put 1:nrow(m), put too big atm
    # browser()
    IDTBPNEW<-IDTBPNEW+1
    subdtNew<-subdtb[subdtb$choice=="TRUE",] # dataframe to replace
    subdtNew$IDTBP<-as.numeric(IDTBPNEW)
    subdtNew$Teaser_rate_years <-m$Teaser_rate_years[i]
    subdtNew$Fee<-m$round_Fee[i]
    subdtNew$Max_LTV<-m$round_Max_LTV[i]
    subdtNew$initial_gross_interest_rate <-m$initial_gross_interest_rate[i]
    subdtNew$BANKID  <-m$BANKID[i]
    subdtNew$choice<-"FALSE"
    subdtN<-rbind(subdtN,subdtNew)
  }
  
  # browser()
  subdtNew<-rbind(subdtb,subdtN)
  
  
  TM2NEW.l<-mlogit.data(subdtNew, choice = "choice", shape = "long", chid.var = "borrowerID", alt.var = "IDTBP")
  
  
  # browser()
  #chose borrowers that chose the same contracts and swap his choice to the new products
  dd<-as.data.frame(TM2NEW.l[TM2NEW.l$choice=="TRUE",])
  a<-summary(as.factor(dd$IDTBP))
  ID_ML<-as.numeric(row.names(as.matrix(a[1]))) # get the numeric value of the most likely occurrence f IDTBP
  
  IDTBPNEW<-max(as.integer(as.vector(unique(TM2$IDTBP))))
  for (i in 1:nrow(m)) {
    IDTBPNEW<-IDTBPNEW+1
    dd[dd$IDTBP==ID_ML,]$Teaser_rate_years[i] <-m$Teaser_rate_years[i]
    dd[dd$IDTBP==ID_ML,]$Fee[i]<-m$round_Fee[i]
    dd[dd$IDTBP==ID_ML,]$Max_LTV[i]<-m$round_Max_LTV[i]
    dd[dd$IDTBP==ID_ML,]$initial_gross_interest_rate[i]<-m$initial_gross_interest_rate[i]
    dd[dd$IDTBP==ID_ML,]$BANKID[i] <-m$BANKID[i]
    dd[dd$IDTBP==ID_ML,]$IDTBP[i] <-IDTBPNEW
  }
  
  dfNEW.l<-dd
  
  return(list(TM2NEW.l,dfNEW.l,nrow(m)))
}
#############################

############################ #to change line in which  R_n <-Matrices[[6]]-0.02
solve_counterfactual_2018<-function(list_a,alpha_d,mean_alpha,TM2,df){
  
  
  timet<-1 # to change when I have more time period
  maxiteration<-1
  
  X_Rnn<-0
  R_n<-10
  iteration<-1    # count the number of iterations
  #browser()
  
  #create variables for the bank that changed its product
  j<-TM2$BANKID[1]
  
  #contraction mapping iteration
  while (max(abs(X_Rnn-  R_n))>0.01 ) {
    
    #get the lender that changed its product
  
    
    # did this for the simulation
    
    if (iteration==1) {
      listaex<-list(list_a[[j]],list_a[[j]],list_a[[j]])
    }else{
      listaex<-list(my.list.expected.profits[[j]],my.list.expected.profits[[j]],my.list.expected.profits[[j]])
    }
    
    
    # lista<-list(my.list.expected.profit.alpha[[1]][[j]],my.list.expected.profit.alpha[[1]][[j]],my.list.expected.profit.alpha[[1]][[j]])
    
    # length(sort(unique(TM2$IDTBP[TM2$BANKID==5&TM2$choice=="TRUE"])))
    
    #do this for all banks
    Matrices<-create_Matrices_2018(j,listaex,alpha_d,mean_alpha)
    
    tilde_LambdaD_jt <-Matrices[[1]]
    tilde_GammaD_jt<-Matrices[[2]]
    
    Lambda_jt_n <-Matrices[[3]]
    Gamma_jt_n<-Matrices[[4]]
    Q_n<-Matrices[[5]]
    
    R_n <-Matrices[[6]]
    # 
    One_M<-  Matrices[[9]]
    #print(j)
    
 
    
    mc_LTV<-unlist(Matrices[[8]], use.names=FALSE)
    fees<-unlist(Matrices[[10]], use.names=FALSE)
    mc_teaser<-unlist(Matrices[[11]], use.names=FALSE)
    mc_fees<-rep(fees[1], length(mc_teaser))
    mc_firm<-rep(j, length(mc_teaser))
    
    mcdt<-data.frame(mc_LTV,mc_fees,mc_firm,mc_teaser)
    Marginal_Costs<-predict(Marginal_Costs_regression,mcdt)
    
    
  #  browser()
    

    X_Rnn<- solve(Lambda_jt_n) %*% ( Gamma_jt_n %*%   R_n + ( tilde_LambdaD_jt - tilde_GammaD_jt  )  %*%   Marginal_Costs- Q_n  %*%  One_M    )
    
    # 
    #update the rate of the products
    R_n_counter<-Matrices[[7]]  
    
    product_data<-unique(df[,c("IDTBP","Fee","Max_LTV","initial_gross_interest_rate","Teaser_rate_years","BANKID")])
    data.rate<-data.frame(as.vector(X_Rnn),sort(product_data$IDTBP[product_data$BANKID==j]))#new rate with their IDTBP
    
    
    # browser()
    #do it just for one counterfactual for now
    for (c in 4) {
      TM20ld<-TM2
      #update the rates
      if (!is.na(R_n_counter[[c]][1]) ) {
        for (idtbp in data.rate[,2]) {
          TM2$initial_gross_interest_rate[TM2$IDTBP==idtbp]<-data.rate[data.rate[,2]==idtbp,1]
          df$initial_gross_interest_rate[df$IDTBP==idtbp]<-data.rate[data.rate[,2]==idtbp,1]
        }
      }
    }
    
 
  
  #update the MS loan size ect (takes some most of the iteration time)
  
  tmp<-create_list_expected_profits_2018(TM2,df) 
  my.list.expected.profits<-tmp[[1]]
  
  
  # browser()
  # print(X_Rnn)
  # print(R_n)
  print( max(abs(X_Rnn-  R_n)))
  #print( mean(abs(X_Rnn-  R_n)))
  #print(iteration)
  iteration<-iteration+1
  
  if (iteration>=maxiteration) {
    return(list(X_Rnn,my.list.expected.profits,df,TM2,0))# 
  }
  
  
  }
  
  return(list(X_Rnn,my.list.expected.profits,df,TM2,1))
}
############################
